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In the context of blackhole perturbation theory, we describe both exact evaluation 
of an asymptotic waveform from a time series recorded at a finite radial location 
and its numerical approximation. From the user's standpoint our technique is easy 
to implement, affords high accuracy, and works for both axial (Regge- Wheeler) and 
polar (Zerilli) sectors. Our focus is on the ease of implementation with publicly 
available numerical tables, either as part of an existing evolution code or a post- 
processing step. Nevertheless, we also present a thorough theoretical discussion of 
^ asymptotic waveform evaluation and radiation boundary conditions, which need not 

Qi be understood by a user of our methods. In particular, we identify (both in the 

^ \ time and frequency domains) analytical asymptotic waveform evaluation kernels. 



and describe their approximation by techniques developed by Alpert, Greengard, and 
Hagstrom. This paper also presents new results on the evaluation of far-field signals 
for the ordinary (acoustic) wave equation. We apply our method to study late-time 
decay tails at null-infinity, "teleportation" of a signal between two finite radial values, 
and luminosities from extreme- mass-ratio binaries. Through numerical simulations 
with the outer boundary as close in as r = 30M, we compute asymptotic waveforms 
with late-time decay {£ = 2 perturbations), and also luminosities from circular 
and eccentric particle-orbits that respectively match frequency domain results to 
relative errors of better than 10~^^ and 10~^. Furthermore, we find that asymptotic 
waveforms are especially prone to contamination by spurious junk radiation. 

PACS numbers: 04. 25. Dm (Numerical relativity), AMS numbers: 41A20 (Approximation 
by rational functions), 44A10 (Laplace transform), 65D20 (Computation of special functions, 
construction of tables), 83-08 (Relativity and gravitational theory. Computational methods), 
83C57 (General relativity. Black holes). 
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I. INTRODUCTION 



Asymptotic-waveform evaluation (AWE) is a long standing challenge in the computa- 
tion of waves. Whether for acoustic, electromagnetic, or gravitational waves the goal is 
to identify the far-field or asymptotic signal radiated to future null infinity J'^ using only 
knowledge of the solution on a truncated, spatially finite, computational domain. For the 
unit-speed ordinary wave equation, the far-field signal is f{u,9,(j)) = limj._j.oo r, 6*, 0), 
where ip{u,r,9,(j)) is a solution to the wave equation written with respect to retarded-time 
u = t — r and spherical polar coordinates {r,6,(j)). Continuing with this example, we view 
computation of roQip{u,roo,0,(f)) as AWE, so long as the evaluation radius r^o can be taken 
arbitrarily large, even if ultimately finite. Indeed, for the ordinary wave equation the asymp- 
totic signal and signal at Too = 10^^ are identical to about double precision machine epsilon 
(see Appendix |A] for the error estimate). The asymptotic signal we compute corresponds 
to observation on the timelike hyper-cylinder r = Toq. In the perturbative gravitational set- 
ting considered here, with r^o = 2M (1 x 10^^) in terms of the Schwarzschild mass M, such 
observations take place in the astrophysical zone pJi3|]. An observation in the astrophysical 
zone is better approximated as taking place at J^"*" rather than future timelike infinity i~^ 
For this reason, and for clarity of exposition, throughout this paper we refer to such 
observations as taking place at J^^. 

Identifying the gravitational wave signal at J^~^ is of both theoretical and practical im- 
portance. Theoretically, in the context of asymptotically flat spacetimes Sachs identified 
the asymptotic metric factors corresponding to f{u, 6, 0) for gravitation, and exploited this 
identification in his discussion of the radiative degrees of freedom for general relativity. 
Using Geroch's calculation framework |6| , Ashtekar and Streubel expanded on the Sachs ap- 
proach in their fundamental analysis |7| of the symplectic structure of radiative modes and 
gravitational flux in general relativity. Several works j8-l3| then investigated the general 



charge integral (where the integration is over a two-surface "cut" of ^+) corresponding to 
the Ashtekar-Streubel flux. The practical importance of AWE stems from the upcoming 
generation of advanced-sensitivity ground-based gravitational wave interferometer detectors 
(i.e., advanced LIGO, advanced Virgo, and KAGRA) [ll-14| and anticipated space-based 



detectors like LISA [15|, ll6||. These instruments are well-modeled as idealized observers 
located at J^^. 

While the problem of AWE for gravitation shares difficulties with its counterparts for the 
ordinary wave and Maxwell equations, for example the slow fall-off of the waves (in our case 
metric perturbations) in powers of 1/r, the gravitational problem is further complicated by 
the backscattering of waves, coordinate (gauge) issues, and non-linearities. For a perturbed 
Schwarzschild blackhole, covariant and gauge invariant approaches exist for the construction 
of "master functions" from the spacetime metric perturbations (see, for example, (l7l-19|). 



In the asymptotic limit these master functions specify the gravitational waveform. Here we 
consider the wave equations which directly govern these master functions. In this pertur- 
bative setting, AWE is precisely the technique (perhaps extrapolation, for example) used to 
compute the master functions at arbitrarily large distances from the central blackhole. This 
paper introduces a new technique based on signal teleportation between two finite radial 
values. 

A straightforward and longstanding approach to AWE in both full general relativity 



20|, |2l| and perturbative settings [22[ has been to record relevant field quantities at a variety 



of radii, perform a numerical fit, and then extrapolate to larger radii. However, the accuracy 
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of this method ultimately relies on an Ansatz for the expected fall-off of the field with larger 



r, as well as recording field values at multiple and preferably large values of r |23[. An 
alternative approach, known as Cauchy characteristic extraction, is to record geometrical 
data from a Cauchy evolution on a world-tube, which is later used as interior boundary data 
for a second characteristic evolution whose coordinates have been compactified to formally 



include J^~^ within the numerical grid |24l-l27l|. When the background coordinates are fixed. 



J^"*" can be directly included within a Cauchy evolution by a geometric prescription using 
hyperboloidal methods 0, 23"34|- Another approach due to Abrahams and Evans shows 



how one may exactly evaluate asymptotic waveforms from gravitational multipoles for general 



relativity linearized about fiat spacetime [35|, |36j. This paper presents a new analytical and 



numerical method to evaluate asymptotic gravitational waveforms from perturbations of 
a non-spinning (Schwarzschild) blackhole. It also presents new results on the asymptotic 
signal evaluation problem for the acoustic (i.e. ordinary) wave equation. Our approach is 
most similar to that of Abrahams and Evans. Essentially, we reformulate their approach 
in a way which subsequently generalizes to blackhole perturbations. However, while our 
approach generalizes the Abrahams-Evans one to a curved background spacetime, we do not 
match their careful discussion of gauge issues. 

For the Einstein equations linearized about Minkowski spacetime in the Lorenz gauge 
the trace-reversed metric perturbation obeys the fiatspace (ordinary d'Alembertian) tensor 
wave equation. Therefore, these perturbations are akin to solutions of either the ordinary 
wave equation or the Maxwell equations; solutions characterized by the sharp Huygen's 
principle and, therefore, which possess secondary lacunae [s^: given trivial initial data and 
an inhomogeneous source which is bounded in space and time, the solution vanishes on the 
intersection of all forward light cones whose vertices sweep over the support of the source. 
The secondary lacunae is a region of spacetime which is "dark" because all waves have 
already passed. Similar statements hold for the homogeneous case with non-trivial initial 
data of compact support. Actually, for the Maxwell case with certain sources, the solution 
may have a quasi- lacunae featuring a late-time static electric field (37| . 

Wave propagation on a curved spacetime is more complicated due to the backscatter- 
ing of waves off of curvature. Even within the relatively simple setting of perturbations 
of Schwarzschild blackholes, backscattering effects are present and the resulting late-time 



"tails" [38|, |39| have been extensively studied both theoretically and numerically. Backscat- 
tering confounds our intuitive sense of "outgoing" and "ingoing" ; one might reasonably take 
the viewpoint that a partially backscattered wave has both outgoing and ingoing pieces. 
Nevertheless, for the linear master equations which describe perturbations of Schwarzschild 
blackholes, there is an unambiguous notion of "outgoing" , provided initial data of compact 
support. Away from the support of the initial data, Laplace transformation of a master 
equation ([T]) yields a homogeneous second-order ODE, which therefore has two linearly in- 
dependent solutions, \l/*^^''(s, r) and '^^~^\s,r), where s is Laplace frequency. Here we have 
suppressed harmonic indices {£, m) and assumed that the area radius r is the independent 
spatial variable. We may assume that \l/(''=^)(s, r) ~ exp(=|=sr=K) as r^r^ — )■ oo, where r* is 
the Regge- Wheeler tortoise coordinate defined below. At a radial location beyond the sup- 
port of the initial data, the frequency-domain solution has the form a(s)\l/^"'"^^(s, r), where 
the details of the initial data are buried in the coefficient a{s). Physically, this notion of 
"outgoing" would perhaps be better characterized as "asymptotically outgoing" . Neverthe- 
less, provided the solution has this form, we can derive at a finite radius both time-domain 
boundary conditions j41-4^ and an AWE procedure. The strategy in both cases is to 
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write down the exact conditions/procedure in the frequency domain, and then accurately 
approximate this exact relationship in a fashion that allows for simple inversion under the 
inverse Laplace transform. For the case of boundary conditions, one approximates the exact 
Dirichlet-to-Neumann map as a rational function (in fact a sum of simple poles) along the 
axis of imaginary Laplace frequency (the inversion contour). The exact time-domain bound- 
ary condition is a history-dependent convolution, which maybe approximated to machine 
precision as a convolution involving a kernel given by a small sum of exponentials. As we 
show, this type of kernel effectively localizes the history dependence. 

Our reformulation of the Abrahams-Evans procedure (and its generalization to curved 
spacetimes) features a similar history-dependent convolution involving a sum-of-exponentials 
time-domain kernel. Section [XTTl demonstrates that, in the (Laplace) frequency domain, an 
AWE kernel is exactly expressible as an "integral over boundary kernels" , thereby allowing us 
to leverage existing codes and knowledge for generating and approximating boundary kernels. 
While the construction of AWE kernels is computationally intensive, this is an offline cost. 
Once the kernel has been calculated, efficient and accurate AWE can be implemented within 
an existing code in a non-intrusive manner. Furthermore, AWE can be effected as a post- 
processing step on existing data recorded at a fixed radial location. Kernels used in this 



paper, as well as others, will be available at [44 



This paper is organized as follows. Section [IT] provides a self-contained guide on using 
boundary and AWE kernels in either existing codes or data post-processing. Towards this 
end. Section 111 CI considers the numerical evolution of late time tails from an approximate 
asymptotic signal which we find to decay at the rate predicted for 1 = 2 perturbations at 
. Section UTTI presents the theoretical underpinnings of both radiation boundary condi- 
tions and waveform teleportation, considering both wave propagation on fiat (Minkowski) 
and Schwarzschild spacetime. For £ = 2,3, 64 perturbations. Section IIV Al considers ac- 
curate signal teleportation to a finite (near-field) radial value. In Section IIVBI we apply 
our method to compute gravitational waveforms and luminosities from extreme mass ratio 
binary systems, finding excellent agreement with frequency domain computations. In these 
studies we have observed that spurious junk radiation is problematic for accurate com- 
putations. Finally, we conclude in Section |V] by discussing open issues, both theoretical and 
practical. 



II. IMPLEMENTATION HOW-TO GUIDE 



Our aim in this section is not to give a derivation of our AWE method. Rather, adopting 
the simplest possible evolution scheme and coordinates, we focus on how AWE is imple- 
mented. By presenting an implementation of our AWE method for a simple scheme, we 
hope to convey the key points to the reader, who will then grasp how to implement the 
method within their own evolution scheme. Since our implementation of AWE relies on cer- 
tain radiation boundary conditions (RBC), themselves essential when working on a spatially 
finite domain, we first describe how to implement these within our simple scheme. Here we 
do not discuss our RBC and AWE methods for different background coordinate systems 
(i.e. Kerr-Schild or hyperboloidal foliations), but return to this issue in Appendix O 

Multipole gravitational perturbations of a Schwarzschild blackhole are described by the 
Regge- Wheeler (axial) and Zerilli (polar) formalisms 45|, |46| . In geometric units the corre- 
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spending "master" wave equations have the form 



dt"^ drl 



+ V'^'^^\r)^! = S, (1) 



where S{t, r) is a possible source and in terms of the blackhole mass M the Regge- Wheeler 
tortoise coordinate is r,, = r + 2M log(^M~^r — 1), which we also denote by x. Until Section 
IIVBI we alway choose 5 = for the source. Expressions for the Zerilli V^{r) and Regge- 
Wheeler V^{r) potentials are given in, for example, Eqs. (2) and (3) of Ref. [4o[. Both 
potentials depend on the orbital angular index i. We have suppressed this index on V{r), 
as well as the orbital and azimuthal indices (£, m) on the mode \l/ and source 5*. 



A. Simple evolution algorithm with RBC 



Most numerical schemes for evolution of ([T]) employ some form of "auxiliary variables", 
for example the variables 11 = —dt"^, $ = dx"^- However, for transparency here we use 
"characteristic variables" : 

^, w = -ii-^, x = -n + $, (2) 

for which the evolution equations are 

d^^if = ^{W + X), dtW = -dxW -V{r{x))'i', S^X = - \/(r(x))^. (3) 

These equations show that W propagates left-to-right {/^) and X propagates right-to-left 
(\). While from the theoretical standpoint we prefer to view fields, for example \l/(t,r), 
as depending spatially on r, from a computationally standpoint we discretize Eqs. ([3]) in 
X. Therefore, suppose the computational domain is the interval [a, b] in tortoise coordinate 
X, corresponding to the interval [ra,rf,] in Schwarzschild radius r. Boundary conditions 
must be specified for VI^ at a; = a and for X at x = b. We discretize ([3]) using first-order 
upwind stencils in space, and the forward Euler method in time. Respectively, let and 
{xk}f=Q denote uniformly spaced temporal and spatial grid points such that Ax = Xk+i — Xk- 
Moreover, let Vk = V{r{xk)) and ~ r(xfc)), with the same notation for the other 
fields. Then the scheme for updating the fields at all Xk from time tn to time tn+i = tn + At 
is given by Algorithm [TJ 
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Algorithm 1 Finite difference scheme for Eqs. 
1: = + {At/2) [W^ + X^] 

2: Wq~^^ = W^~^^ > Boundary condition at a 

3: X^+^ = Xo" + (At/ Ax) (Xr - Xo") - At^o^o 

4: for A; = 1 to K - 1 do 

5: ^^+1 = ^l+ {At/2) [W^ + X^] 

6: W^^' = - {At/ Ax) {W^ - W^_,) - AtV^n 

7: ' = X]: + {At /Ax) {X^^, - X^) - AtV^n 

8: end for 

9: ^^+1 = ^l + {At/2) {W^ + X)^) 
10: Wp-^ = W^- {At/Ax){Wp^ - W^_^) - AtVK^]^ 

11: Xp'^ = X^+^ > Boundary condition at b 



To complete the scheme, we must specify both ~ W{tn,ra) and X" ~ X(t„,r;,) as 
functions of (discrete) time. So long as a ^ 0, the inner boundary value of the potential 
V{r) near r = is zero to machine precision, and the Sommerfeld boundary condition 
= 0, Vn is highly accurate. The RBC at x = 6 is determined by a Laplace convolution, 

X{t, n) = IM f E{t - t', r,)^{t', n)dt\ (4) 
rb Jo 

where f{r) = 1 — 2M/r and the boundary time-domain kernel is 

m,n) = t'^^E,(t,n), z,(t,n) = e.J^) ^ (5) 

g=l ^ ^ 

The parameters {{'yq{pb) , (3q{Pb))}q=i depend on the rescaled boundary radius ph = {2M)~^rb 
(as well as the orbital index i which is suppressed) and they are listed in numerical tables, 
such as those given in the Appendix [pfl Some of the parameters {{jq, Pq)}q=i are complex, 
but the kernel S(t, r?,) is real. We stress that, insofar as implementation of the convolution 
is concerned, the origin of these numbers is unimportant. Defining, for example, '^b{t) = 
\E'(t, Tfe), we write (jl]) as Xb{t) = r^^/(rb)(S(-, r^) * ^6)(t)- With a similar notation, the 
constituent convolution {Eg * \E'5) = {Eq{-,rb) * '^b){t) obeys an ODE at the boundary, 

|(H, * VP,) = [{2M)-'(3,{E, * VI/,) + V]/,] . (6) 

These d ODE can be integrated along side the system ([3]). Indeed, we define {Eg * vE^b)" ~ 
{Eg * v]>,)(t^) and complete our scheme as follows. 



These parameters can be redefined through division by 2Af , thereby removing 2M factors in the formulas 
which follow. Indeed, such a redefinition would be made in an actual code. Nevertheless, we retain the 
original {{'Jq{pb) : f3q{pb))}q=i parameters with 2AI factors, in order to ensure that the parameters here 
correspond precisely to those listed in our tables. 
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Algorithm 2 Finite difference scheme for Eqs. with RBC 

1: for g = 1 to d do > First update constituent convolutions 

2: (S, * ^5)"+' = (Hg * ^bT + [(2M)-i/3,(Sg * + 

3: end for 

5: W^+^ = 

6: Run Algorithm [T] with PV^^"*"^ and X^~^^ given by above values. 



B. Evaluation of the asymptotic waveform 

Our approach to AWE is similar to the described implementation of radiation boundary 
conditions. We introduce new parameters {{'jgipb, Poo), f^gipb, Poo))}q=i and a new kernel 

^Eu . ^ lg{pb,Poc) ^Eu \ \ (PliPlnP^\ 

^ (t,rb,roo) = 2^ — (t,p6,Poo), ^Jt,Pb,Poo) = exp 1^ — j. (7) 

The parameters {(7^(pb, Poo), f^qiPb, Poo)}^^ 1 also depend on i, but this dependence has been 
suppressed. The E here stands for "evaluation" and differentiates this kernel from the RBC 
one. This evaluation kernel enacts teleportation (the term is defined in Section IIIII) of the 
waveform from the boundary radius = 2Mpb to the evaluation radius Too = 2Mpoo ^ 
(from b to Xoo in tortoise coordinate). As discussed below, the choice Too = 00 is formally 
possible (see Sec. |V]); however, in this paper is an arbitrarily large, albeit finite, radius. 
To ensure the signals recovered at Too and J^"*" are identical to about machine precision, 
we choose Too = 2M(1 x 10^^) for double precision simulations, and would choose = 
2M(1 X 10^°) for quadruple precision simulations. We then approximate the asymptotic 
waveform as 



^oo(t) :=^(t+(xoo-&),roo) =^ / S^(t-t',n,roo)^(t', 

^0 



rb)dt' + ^{t,rb). (8) 



The offset by \E'(t, rj,) in this formula stems from a technicality explained in Section [III A[ 
As before, this formula can be implemented through integration of ODE at the boundary, 
only now these ODE are not coupled to the numerical evolution. With (S^ * '^b)"' — 
(S^(-,rf,) * \t'fe)(t„) and "if^o — ^oo(^n), the algorithm using forward Euler is as follows. 

Algorithm 3 Evaluation of asymptotic waveform, placed after field update 



for q = 1 to do 

(Sf * = (Sf * + At[(2M)-i/3f (Sf * ^,)" + 

end for 



The following condition would yield particularly efficient AWE: 

= d, (iqiPb, Poo) = fiqipb), Vg, (preferred, but perhaps not possible). (9) 
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FIG. 1. QuASiNORMAL RINGING AND DECAY TAILS. Each dashed curve corresponds to power law 
decay, with the indicated rate determined by a least squares fit of the field over the time window 
[500, 700]M. The time shift for the asymptotic waveform has not been included. 

Indeed, integration of the same ODE (E]) at the boundary would then determine both the 
RBC and AWE. In this case, steps 1 through 3 of Algorithm |3] have already been carried 
out in Algorithm |2l However, the assumption in (|9]) may not always be possible, and even 
when possible appears to yield less accurate teleportation kernels. We have constructed 
^ = 2 AWE kernels which satisfy (|9]) (an example is given in the Appendix |D]); however, 
relative to our best kernels, they indeed yield less accuracy. Moreover, we have been unable 
to achieve (jH]) for £ = 64 teleportation kernels. Therefore, in this paper we will not assume 

(ED. 

Remark. If ^ is itself complex (as is the case in applications), then round-off issues will 
lead to mixing of the real and imaginary parts in the simple algorithms above. In this case, 
we advocate splitting the complex exponentials which make up S and into manifestly real 
expressions involving sine and cosine terms. Such splitting amounts to extra bookkeeping, 
but hardly complicates the above treatment. 



C. Numerical experiment: quasinormal ringing and decay tails 

Section HV] documents the results of several numerical experiments which validate and 
test our methods. This subsection also describes a numerical experiment, although here 
with the goal of providing further assistance toward implementation. An interested reader 
might first repeat the experiment described below. 

We consider the £ = 2 Regge- Wheeler equation and Gaussian initial data 

^^^-[2-./(2M)P^ ^^'^^_^^-[2-.K2M)]\ n = <|.(0,r(x)). (10) 



For this experiment the computational x-domain is [a, 6], where a = — 200M and b = 
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5000 10000 15000 20000 25000 

t/M 

FIG. 2. Tail decay rate for the teleported signal. The rate p for ^oo{t) oc has been 
computed using logarithmic difference quotients based on d\^t\n.\^ oo{t)\. 



30M + 2Mlog(14) corresponds to pb = n/{2M) = 15. We evolve the data until time 
t = 600M, with a Sommerfeld boundary condition at x = a and the convolution boundary 
condition (jlj) at x = 6. Table IIIII in Appendix [D] lists the 26 pole RBC kernel S used for this 
convolution. Instead of Algorithm [T] we use a multidomain Chebyshev collocation method 
with classical Runge Kutta 4 as the timestepper. 

During the simulation we record as a time series both the field \l/()f:,20M) and \l/(t,rf,). 
Through the convolution ([8]) determined by Table IVl in Appendix [D| the field \l/(t,rf,) is 
teleported from = 30M to Too = 2M(1 x 10^^) providing a time series which approximates 
the asymptotic waveform \l/oo(t). In absolute value \l/(t,20M) (solid blue line) and \l/oo(^) 
(solid black line) are depicted in Fig. [TJ The time shift for the teleported waveform has not 
been included, and we have chosen to record \l/(t, 20M) at r = 20M < r?, only to ensure 
that the time series in the plot do not lie on top of each other at early times. These series 
exhibit the phenomena of quasinormal ringing and late time decay tails. Each dashed curve 
in the figure corresponds to power law decay, with the indicated rate determined by a least 
squares fit based on the numerical decay of the field over dashed curve's time window. The 
decay rates p = —7 and p = —4 are respectively the theoretical predictions (sol . 47 . for 
a finite radius and 

Figure [2] shows the decay rate computed from p = 9intln|\l/oo(t)| for the teleported signal 
with the evolution carried out to t = 25000M. As a post-processing procedure, generation 
of this figure from existing data would take a few seconds. The decay for the teleported 
signal asymptotically approaches p = — 4 at later times. The teleported pulse corresponds 
to a time series recorded along the wordline (tobs + — b,r^) by an observer safely within 
the astrophysical zone. Here x^o — b is the time shift, and the astrophysical zone is defined 
as the region where tobs ^ with tobs the time elapsed after the pulse's leading edge 
passes our fictitious observer [l|-l3|. Since an observation in the astrophysical zone is well 
approximated as taking place at J^"*", the observed —4 decay rate is expected. For very 
late times the decay rate should settle towards —7, although extended precision might be 
necessary to capture the transition. 



III. THEORETICAL DISCUSSION 



To fix ideas and motivate the new method, the next subsection describes AWE for 
fiatspace multipole solutions of the ordinary 3+1 wave equation, formally the M = case of 
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Eq. (11]). Formulas derived in the next subsection motivate similar ones given for blackhole 
perturbations in subsection IIII B[ 



A. Flatspace waves 



This subsection describes (i) outgoing multipole solutions to the flatspace radial wave 
equation, (ii) the exact RBC obeyed by outgoing multipoles, and (iii) the relationship be- 
tween this RBC and AWE for outgoing multipoles. Throughout this subsection, we of- 
ten choose £ = 2 as a representative example, but similar (obvious) results hold for any 
multipole-order i. 



1. Structure of outgoing and ingoing flatspace multipoles 

General outgoing (e = 1) and ingoing (e = —1) order-£ multipole solutions of the 3+1 
wave equation {—d'^ + d1 + dy+d1)il) = Q have the form (see, for example, Refs. 47, [HI ]) 



^(t,a:,y,z) = -vl/(^)(t,r)y,^(^,0), ^f{t,T) = J^'-Cekf^'-'Ht-er), c,, = 

(11) 

where we have suppressed the azimuthal index m on the "mode" ^5'^(t,r). In Eq. 
f^^\u) is the pth. derivative of an underlying function f{u) of retarded time u = t — r, and 
similarly for f^^\v), where v = t + r is advanced time. In (ITTl) we view (x,y,z) as place 
holders for (rsin^cos(/>, rsin^sin0, r cos6'). For both e = ±1 cases the mode obeys the 
flatspace radial wave equation 

and, specializing to the representative example, the outgoing quadrupole is 

¥^\t, r) = fit - r) + -fit - r) + ^/{t - r). (13) 

We are interested in the (e = 1) outgoing case, and so now write \l/^(t,r) to mean \I'^^\t,r). 

Given a fixed radius which specifies the outer boundary, consider the following assump- 
tion on the initial data (one of compact support): 

^'<;(0,r) = = ((9t*<?)(0,r), r>n-5, for any small 5 > 0. (14) 

Provided that f|T4|) holds, in the region r > rt, Laplace transformation of an outgoing mode 
\&£(t, r) yields 

e 

$,(., r) = a{s)s'e-We{sr), W,{z) = ^, (15) 

k=0 
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where z = sr and 

pOQ poo 

a{s) = e'^ / e-^V(t - r)dt = / e-'''f{u)du. (16) 
Jo J-n 

Notice that a{s) is indeed independent of r. For the quadrupole case, we have 

$2(s, r) = a{s)s\~'''W2{sr), W2{z) = 1 + ^ + ^ (17) 

[compare this expression with Eq. f|T3|) ]. 

To obtain ( 1151] and fll6l) . we have used f[T^ as follows. First consider the Laplace transform 
of /*^^~'^)(t — r) which appears in flTTl) . Repeated integration by parts generates t = 
boundary terms of the form y(^~'^~p)(— r) for 1 < p < i — k. All such terms vanish, as can 
be shown by the following identity: 

fit -r) = (-l)'(lyy (dt + d^)Y^e{t, r). (18) 

Since the initial data and dt^ ^ vanishes on an open spatial neighborhood of the spatial 
point with coordinate r > r^, in fact vanishes in an open spacetime neighborhood of 
the spacetime point with coordinates (0,r). Therefore, all derivatives of '^i vanish in the 
same neighborhood, which implies r) = 0. This implication stems from repeated 

differentiation of (ITSl) with t = enforced afterward. 

The previous argument establishes that 

/•oo /"OO 

/ e-''f^^-''\t-r)dt = s^-^e-''' e-'''f{u)du. (19) 

Jo J~r 

The lower limit — r of integration can now be replaced with — r^. Indeed, f|T^ implies that 
\l/(t, r) = for < t < 5 + (r — Tf,) and r > r^. Using f|T8|) . we conclude that f{u) = for 
—r<u< —Tb- Therefore, the right-hand side of the last equation is a(s)s^~^e~^''". 



2. Radiation boundary conditions for flatspace multipoles 

We continue to derive expressions for r > Vb where the assumption (fT^ of compact 
support holds. The explicit expression ( ITSl) for \l/£(s,r) determines an exact frequency- 
domain boundary condition 

s^i{s, r) + dr^iis, r) = ^Q^is, r)$^(s, r), (20) 

where the frequency- domain radiation kernel f2^(s,r) defines the Sommerfeld residual. In- 
deed, the operator on the left-hand of (l20l) corresponds to the Sommerfeld operator dt + dr 
in the time-domain. If £ = 0, then n^(s,r) = 0; otherwise a simple computation based on 
Eqs. (I15f20p shows that the frequency-domain kernel is given by (with the prime indicating 
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differentiation in argument): 

where {bi^k : I < k < i} are the roots of We{z), all of which are simpleU For the quadrupole 
case 

S.(..r) = ^±^ + ^^. -. = -|±if. (23) 

s — z+/r s — Z-jr 1 1 

where z^ = 62,1 and Z- = 62,2 solve W2{z±) = 0. 

The time-domain RBC is the inverse Laplace transform of (!20l) . i.e. the Laplace convo- 



lution 49-53 



dt^e + dr^e = - f m - t', r)^i{t', r)dt\ r) = V ^ exp f ^) 

r Jo r \ r J 



(24) 



Subject to our assumption (fT^ of compact support, the outgoing multipole [e = 1 in Eq. (fTT!) ] 
obeys Eq. ( !24l) exactly, as can also be shown via direct calculation using repeated integration 
by parts; see Appendix [Bl If, on the other hand, assumption (HM does not hold, then Eq. (12^ 
is violated, but only by terms which decay exponentially fast in t; again, see Appendix |Bl 



3. Asymptotic waveform evaluation and teleportation for flatspace multipoles 

For a generic outgoing solution, it is possible to recover the profile function f{t — r) and 
asymptotic waveform 

^^(t,r) ~ /W(i-0, r^oo, (25) 

via data recorded solely at a finite and fixed radial location, again taken as r = r?,. Let us 
consider the i = 2 concrete example. Generalization to higher i is straightforward. 

Equation ([T^ suggests that we solve the ODE initial value problem 

y" + -y' + ^y = ^2(t,r), |/(0) = = y'iO), (i = 2 problem) (26) 

in which case f{u) = y{u + r). 

In a pioneering series of papers jsS, 3^, Abrahams and Evans showed how the above pro- 



cedure carries over to the theory of gravitational multipoles for general relativity linearized 
about flat spacetime. We now re-examine the basic idea behind Abrahams-Evans AWE from 
the standpoint of Laplace convolution, and will consider two kernels: one for evaluation 
of the underlying function f{u) and another $^ more suited for evaluation of the waveform 
f^^\u). Our implementations have mostly relied on the kernel. 

^ The last equality also follows from the identity 

Wi{z) = J—e'K,+^,2{z), (22) 

V TT 

showing that the h^^k are also the roots of the half-integer MacDonald function if £+1/2 (2:), which are 



simple and lie in the left-half plane 55l. l56l|. The appearance of Ki^i/2{z) may have been anticipated; 
indeed, the modified Bessel equation arises when finding separable solutions to the Laplace transformed 
flatspace radial wave equation ([T2|) . 
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Continuing with the 1 = 2 example, we introduce a frequency- domain 'profile evaluation 
kernel 02(s,r) tailored to satisfy 



62(s,r)\E'2(s,r) = a(s)e~ 



vis)- 



(27) 



That is, the product of 62(3, r) and "^2(3, r) is e f{t — r)dt [cf. Eq. (fT6l) ]. Comparison 
with ( ITTll immediately shows that 



02(3, r) 



1 



ir 

7! 



1 



s^W2{sr 

The corresponding time-domain profile evaluation kernel is 

ir 



(28) 



e2(t,r) 



exp 



r 



exp 



r 



(29) 



and ?/(t) = (02(",'") * ^2(';'"))(^) solves the Abrahams- Evans initial value problem 
Essentially the same arguments show that the order-£ profile evaluation kernel is 



1 



s^WAsr) 



(30) 



Despite appearances, the kernel is regular at s = 0. For example, 1/W2{sr) ~ (sr)^/3, and 
in general l/WAsr) ~ {srY/cge, as s — 0. 

Direct evaluation of the asymptotic waveform is also possible. Teleportation by a positive 
shift r2 — Ti means conversion of \E'(t, ri) to \E'(t + (r2 — ri),r2), and it might correspond 
to a small finite shift r2 — ri. However, when r2 is suitably large, we write r^o for r2 
and view teleportation as an AWE procedure (in which case typically ri = r^, and it is 
the boundary waveform \E'(t,rf,) which is teleported). Teleportation is accomplished with a 
frequency- domain teleportation kernel 



^e(s,ri,r2) = -1 + 



We{sr2) 
WAsrA 



rigged to satisfy 



r2) 



<!>eis,ri,r2)'^i{s,ri) + '^As,ri). 



(31) 



(32) 



We have included the —1 factor in fl^ to ensure that $^(s, ri, r2) has a well-defined inverse 
Laplace transform $^(t, ri,r2). In the time domain we recover the desired property 



(^"2 -ri),r2) = ($£(-, ri,r2) * ^^(■,ri))(t) + ^£(t,ri). 



(33) 



Adjusting for the (r2 — ri) time delay, this formula allows for conversion of the signal at 
ri to the signal at r2. Since r2 < 00, this method can also be used for evaluation of the 
asymptotic waveform f^^\u). We refer to the r2 = 00 case <l'^(s,ri,oo) as the frequency- 
domain waveform evaluation kernel. 

The relationship between RBC and AWE/teleportation kernels is a key insight of this 
paper, and the one which is exploited to numerically construct AWE/teleportation kernels. 
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For example, the profile evaluation kernel can be written as 



exp 



V 



dr] 



(34) 



1/We{sr) 



That the underbraced quantity is indeed 1/Wi{sr) follows easily from the identity ?7^^fi^(s, rj) ■ 
d^logWi{sri), that is essentially the definition f l2T|) . The integration in flM|) can of course be 
carried out, recovering fl28|) for the £ = 2 case; however, when considering similar expressions 
for blackhole perturbations at least some of the integration will be performed by numerical 
quadrature. Similarly, one can express the teleportation kernel through 



^^(■s,ri,r2) = -1 + exp 



r2 



ri 



dr] 



(35) 



Weisr2)/Wi{sri) 



In Section IIIIBI we introduce the analogous kernels for waveform teleportation in the 
Regge- Wheeler and Zerilli formalisms. As mentioned, we have mostly used the kernels. 



4- Efficiency and storage 



Here we comment on RBC and AWE for the ordinary 3+1 wave equation from the 
standpoint of efficiency and storage as £ — )■ oo, both summarizing known results for RBC 
(49I ] and considering these issues for AWE. Let A represent a characteristic wavelength, say 
determined by the initial data or inputted boundary conditions. For numerical evolution 
to a fixed final time T, an implementation of the exact fiatspace RBC f l2^ (with a kernel 
comprised of i exponentials) has the following work and storage requirements: 

WorkexactRBC = 0(A"^£T), Storagee^actRBC = 0{i). (36) 

These scalings are deduced from the cost of integrating i ODE of the form (jH]) with approx- 
imately X~^T timesteps. 

A spatially and temporally resolved numerical integration (with arbitrary boundary con- 
ditions) of Eq. (fT2!) on a radial domain of fixed size corresponds to the following work and 
storage scalings: 0(A~^T) and 0(A~^). Indeed, a resolved spatial discretization of Eq. ( !T2|) 
yields a coupled system of approximately A~^ ODE. As more spatial/temporal resolution 
is typically required for large i solutions, it is reasonable to view A~^ ~ i, in which case 
the scalings for the interior solver are comparable to f l36p . However, implementation of the 
exact RBC is still preferable to choosing the computational domain so large that the outer 
boundary is casually disconnected from the wordline of an interior "detector" . Spatial dis- 
cretization after such domain enlargement yields A~^T coupled ODE, whence 0(A~^T^) and 
O(A^^T) for the work and storage. 



Kernel compression yields a more efficient implementation of RBC. As proven in Ref. (49|, 
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the kernel fle{s,r) admits a rational approximatiorH 



n=l 



E 



s,r) -^i[s,r) 



(37) 



where £ is a prescribed tolerance and the number of approximating poles scales like [49 

d = 0{logu \og{l/e) + log^ u + y-^ \og^{l/e)) (38) 

asz/ = £+ l/2— T-oo and e — t- 0+. The frequency domain bound in (|3711 implies a long-time 
bound on the relative convolution error in the time-domain, see Appendix \M Since d grows 
sublinearly in £ and the approximation S£(s,r) [likewise its inverse Laplace transform 
^^(i:, r)] is called a compressed kernel. An implementation of Laplace convolution RBC based 
on compressed kernels 'Ei{t,r) scales like 



(39) 



WorkcompressedRBC — 0(A ^rfT), Storage^ompressedRBC — 0{d), 

with clear performance in the large-£ limit. 

The proof of relies on the large-£ asymptotics j49|, |55| of the roots {bi^k : k = 1, . . . ,i} 
of Kf^^ /^iz). Precisely, as ^ — )■ oo the scaled roots be^k/{^ + 1/2) accumulate on a curve C 



given by [49|, |55 



AtanhAiiVAcothA- A2, Ag[0,Ao], tanhAo = l/Ao. (40) 



Since the pole locations appearing in both the exact flatspace RBC and AWE kernels are 
{be,k/i^ '■ k = !,...,£}, we conjecture that an implementation of AWE based on kernel 
compression formally satisfies the scalings ( 139|1 . However, we are unsure if these scalings 
hold in practice. 

As a nascent investigation, we consider compressed kernels for £ = 64 flatspace RBC and 
teleportation. Figure |3] plots scaledpole locations for a 20-pole compressed kernel S64(s, 15) 
which approximates flgAis, 15) and for 20, 28, and 36-pole versions of a compressed kernel 
S^(s, 15,240) which approximates <l'64('5, 15, 240). Here we have scaled all pole locations 
by a factor r/(£ + 1/2) = 15/64.5 in order to plot them relative to the curve C, on which 
the actual scaled zeros 664,fc/64.5 lie (at least to the eye). Figure [3] shows that, compared 
with poles for compressed teleportation kernels, the poles for the compressed RBC kernel 
lie much closer to C. Nevertheless, for both compressed RBC and teleportation kernels as 
the number of approximating poles increases (corresponding to a smaller tolerance e), more 
of the approximating poles "lock on" to C. This behavior is evident in the right blow-up 
plot, where for 20, 28, and 32-pole compressed teleportation kernels, we respectively find 
(circles), 1 (diamond), and 3 (squares) "locked-on" poles. Moreover, at least to the eye, 
these correspond to "locked-on" poles (crosses) for the compressed RBC kernel. Sec. IIIIB4I 
briefly discusses these issues for the gravitational case. 



^ The -ji^n and (3i^n appearing in the approximate (frequency-domain) flatspace kernel p7|) are different 
than the similar parameters appearing in the approximate (time-domain) blackhole kernel ([5|). Here je.n 
and /3i,n do not depend on r, whereas the parameters in ([5]) do depend on the (rescaled) radius. We use 
similar notations for the flatspace and blackhole cases, hoping this practice does not cause confusion. 
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FIG. 3. Scaled pole locations for compressed kernels. Here TLP means teleportation, 
and the curve C is described by the parameterization z(A) given in (j40p . In the blow-up plot there is 
a cross-diamond-square coalescence on C (near its right end). See the text for further explanation. 



B. Blackhole perturbations 



We consider the following rescaled versions of the generic master equation ([T]) (retaining 
the same stem letter V for the potentials): 



dr'^ dpi 
here in terms of rescaled coordinates 



i ^ yRW,Z 



(p)^^ = 0, 



p = r/{2M), T = t/{2M), = p + log(p-l) 



(41) 



(42) 



Expressions for the Regge- Wheeler Vf^^(p) and Zerilli V^{p) potentials are given in Eqs. (3) 
and (4) of Ref. 43 1 (expressions in terms of r rather than p are given in (4o|) . The formulas 
we present here hold for both formalisms, and have been drawn from Refs. 4l|-|43| . As before 
in our analysis of the flatspace radial wave equation, here we also suppress the azimuthal 
index m on 



1. Structure of outgoing solutions 

With cr = 2Ms the rescaled Laplace frequency, formal Laplace transformation of ([1]) 
yields 

-^ + Vr''ip)^i + ^'^i = 0, (43) 
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FIG. 4. Profiles for an p = 15, £ = 64 Regge- Wheeler frequency domain RBC kernel. 

Outgoing solutions of the last equation can be expressed as an asymptotic serie^ about 
p = oo, 

^,{a,p) = a{a)a'e-''''*Wi{ap,a), W,{ap,a) ~ 1 (44) 

[cf. Eq. (fT5|) for a flatspace multipole]. Notice that We{ap, a) = Wi{sr, 2Ms), and so formally 
the flatspace expression Wi{sr) = Wi{sr,0). 



2. Radiation boundary conditions 



We again assume initial data of compact support, namely that \E'^(0, p) = = ((9t-\E'^)(0, p) 
for p > pb — 6, where pb specifies the outer boundary. We now work with any p > pf,, in terms 
of which exact radiation conditions satisfied by a generic asymptotically outgoing multipole 
have the following frequency- domain and time-domain forms j41-43|: 



a^e + dp^^e = 1 |^1 - Q^^^ ^ 9,^^ + dp^^e = 1 |^1 - 1^ * (45) 



where the frequency domain radiation kernel is 



ue{a,p) = op 



Wi{(jp, a) 



(46) 



with the prime denoting differentiation in the first argument. 

Refs. 4l|, |43| have argued that the kernel ^^(cr, p) has the following "sum of poles" rep- 



^ The coefficients de^k{cr) defining the asymptotic series We{z, a) ^ J^'kLa ^i-k{o')z are respectively defined 
by three-term and five-term recursion relations in the Regge- Wheeler and Zerilli formalisms 4l|-|43| . 
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FIG. 5. Profiles for £ = 2 Regge- Wheeler frequency domain RBC kernels. As p 
increases these profiles shrink towards the origin (the same phenomena occurs for the flatspace 
RBC kernels). 



resent at ion 



E 



N, 



a 



(P) 



k=l 



a 



Mx; p) 



dx, 



(47) 



where fe{XiP) = p) and the cr^^k{p) are simple roots of Wi{ap,a) (analogous to 

the roots bi^k/r of Wi{sr) in the flatspace case) . At least for p > 15, the integer Ni = i oi 



+ 1, when i is respectively even or odd [41|, |43|. The orig in of the extra root, relative to 



the flatspace case, in the odd-£ case is discussed in Ref. [43|. 

Insofar as numerical implementation is concerned, a key requirement is the ability to 
evaluate the profiles ReQi{iy,p) and lmui{iy,p) for y E M.. These evaluations are along 
the imaginary cr-axis, typically the inversion contour for the inverse Laplace transform. 
Accurate methods for such evaluation have been described in 41 . In fact, these methods 



are not based on the sum-of-poles representation P7|) . but this issue is of no concern here. 
For example, the profiles for an p = 15, ^ = 64 Regge- Wheeler kernel are shown in Fig. HI 
As p increases these profiles "shrink" towards the origin (as do the corresponding flatspace 
profiles), and this phenomena is documented in Fig. [51 With the ability to numerically 
generate the profiles ReQi{iy,p) and Imui^iy, p), we are then able to construct approximate 
kernels via Alpert-Greengard-Hagstrom (AGH) compression [49]]. Here a compressed kernel 
is a sum of simple poles. 



6(^,p) = Yl 



=1 



7£,g(p) 

- Mp) 



0Je{a,p), Re/3i^g{p) < 0,Vg 



(48) 



^ This terminology is suggestive only, since in complex analysis poles are isolated singularities. Therefore, 
the integral term a}™*(f7, p) appearing in (|47|) is not, strictly speaking, a "continuous distribution of poles." 
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FIG. 6. y-GRiD TO EVALUATE THE PROFILES ReD2(iyj,p) AND ImtD2 (ij/j , /)) , Vj . This grid has 21 
points, 4 adaptive levels, and 9 points at the bottom level. 

where the approximation ^i{a,p) satisfies 



with e a prescribed tolerance. The number d clearly depends on e and £, and the numbers 
f3e^q and ■je^q depend both on the boundary radius p (as indicated) and on £ (the dependence 
on which we have restored here). The modifier compressed in the description of C,e{<J,p) is 
apt. Indeed, as described in Sec. IIII A4t for the ordinary wave equation the exact frequency 
domain kernel admits a similar rational approximation with d scaling as ( l38|) . Similar scaling 
has been observed empirically for approximations ^^((j, p) of blackhole kernels ueicr, p) j4l| . 



Algorithm H] summarizes our implementation of AGH compression (see Ref. j42| for a 
complete description). Let us further comment on Algorithm HI with numbers appropriate 
for i = 2. Typically ?/max = 300/p for step 1. For step 2 we have typically chosen 10 to 
20 adaptive levels centered around the origin with 65 points at the bottom level, about 10^ 
grid points yj in all. Fig. E] depicts an example y-giid. For step 3 the evaluation at each yj 
requires ODE integration in the complex plane with upwards of 10^" fioating point operations 
in double precision (more in quad precision). Step 5 is a confirmation step meant to verify 
(149|) . Ideally, this confirmation takes place with a much larger y-window than [— ?/max, 1/max], 
and on a different (dense and uniform) y-grid. This step involves further ODE integration 
and is therefore as or more expensive than step 3. 



Algorithm 4 Steps for compressing an RBC kernel. 

INPUT: i, p = Pbi d (orbital index, dimensionless boundary radius, desired number of poles) 
OUTPUT: {Pi,q{p),li,q{p) : g = 1, . . . , c?} (compressed kernel) 

1: Choose an approximation window [— ?/max, 1/max] on the a = iy imaginary axis. 

2: Partition [— ?/max, I/max] to form a y-grid, typically with mesh refinement at the origin. 

3: Numerically evaluate the profiles Rea}2(i?/j, p) and Ima}2(i?/j, p) on the y-giid. 

4: Compute the numbers {/3^,g(p), 7£,g(p) : g = 1, . . . , c?} by AGH compression. 

5: Using {f3i^q{p) , je^q^p) : q = 1, . . . , d}, verify (H^ . If not verified, repeat with d -i^ d + 1. 



From the standpoint of implementation, the representation fj^8|) is crucial, since it implies 
that the time-domain convolution can be approximately evaluated via integration of ODE 
at the boundary. For a typical explicit ODE scheme and a sufficiently small time-step, 
integration of these ODE is numerically stable since the relevant poles in ( l48l) lie in the 
left-half plane. Let ^^(r, p) be the inverse Laplace transform of ^^(cr, p) with respect to a. 



<^£{(^, P) - ^M, P) <£^p£(o-,p)| 



(T G 1. 



(49) 
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Then the approximate time-domain kernel S£(t,r) = (l/2M)^^(t/(2M),r/(2M)) appearing 
in is the inverse Laplace transform (with respect to s) of 

%is,r) = |:^H,,(.,r), g,,(.,r) = ^ _ (2m1-i/3,,(p) ' ^^^^ 
where, unlike in Section [Til here ^-dependence has not been suppressed. 



3. Asymptotic waveform evaluation and teleportation 



Similar to before, we introduce two kernels: (i) one 6i for evaluation of an underlying 
profile, and (ii) another (f)£ for AWE/teleportation of the waveform. The first type of kernel 
is defined by [cf. Eq. ([3] 



a 



V 



drj 



(51) 



and satisfies 6i{a, p)'$e{a, p) = a((T) exp(— crp*), as can be seen directly from Eq. (14^ . We 
can not analytically perform the integration here. However, the pole part of the kernel can 
be exactly integrated to remove the singularity. Indeed, we find 



p) 



exp 




exp 




exp 


p ^ . 





(52) 



Teleportation (pi — )■ P2 < oo) is defined through the kernel [cf. Eq. ([35 



0£(cr, P1,P2) 



-1 + exp 



pi 



V 



-dr] 



(53) 



Wt (o-p2 ,o-)/Wt (crpi ,a) 



Adjusting for the (p2 — pi) time delay, this kernel allows for conversion of the signal at pi 
to the signal at p2, as it satisfies 



e-(/'2 Pl)^^{a, P2) = 4>i{a, pi, p2)$^(a, pi) + $^(a, pi). 
In the time domain we therefore recover the desired property 

^,(r + (p; - p*), P2) = {M-,Pu P2) * ^£(-, Pi))(r) + ^,(r, pi). 



(54) 



(55) 



Exact evaluation of the asymptotic waveform corresponds to p2 = 00. 

Let us describe how we numerically approximate (j)e{a, pi, P2) as a pole sum ^^(cr, pi, P2). 
A simple version of the procedure is to follow the steps listed in Algorithm [U replacing step 
3 with evaluations of the profiles Re0^(iyj, pi, P2) and lm(j)i{iyj, pi, P2). To generate these 
profiles, we use ([53]) . and for each yj evaluation point perform the r] integration using a com- 
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posite Gauss-Kronrod rule. Assuming [A,B] is one subinterval of [pi, P2] in the composite 
rule, the details are as follows. Let r] = AB/[{B — A)q + A] for q G [0, 1], so that 

^MMt,,^"f ^M^ (56) 

where {qp, u'p)^f are the 15 nodes and weights for the Gauss-Kronrod rule relative to [0, 1], 
and by the above equation we mean the same rule is separately applied to the real and 
imaginary parts of uji{iyj,ri). This numerical integration is accurate, because, in fact, for 
each grid-point Uj integration of both Rea}^(i?/j, r^) and Imujiiiyj^rf) always involves terms of 
the same sign. With Nc denoting the number composite subintervals, profiles for the RBC 
kernels on the y-grid must be computed Nqk-Nq times. If p2 = 00 the last composite interval 
might be handled through a semi-infinite quadrature. Approximation of 6£{a,p) follows a 
procedure similar (although more complicated) to the one outlined above for ^^(cr, pi, ^2)- 

Unfortunately, this simple procedure becomes too costly when p2 ^ pi- The problem 
is two-fold. First, Nc must be chosen large. Second, and more serious, the approximation 
window [— i/max, 2/max] IS fixcd by the profiles for pi (the "widest" profiles in the integration). 
However, since, as seen from Fig. [5l the RBC profiles shrink as p increases, the y-grid needs 
many adaptive levels to resolve the contribution to the ?7-integration from the profiles at 
and near p2- The y-grid then must have a large number of points, on which Nqk ■ Nc 
complex function evaluations are made (with each such evaluation costing ~ 10^ floating 
point operations). To bypass this issue, we follow another rather complicated procedure, 
whereby the interval [pi,P2] is first broken up into Af chunks [pa, Pa+i] which are typically 
decades like [10", 10°"'"^] if p2 is very large. Each chunk has its own approximation win- 
dow [— 2/max5 Z/max] (relatively small) |/"-grid, and we choose these to conform with how 
shrunken the profiles are for p ~ po,. Next, using the relatively simple procedure described 
in the last paragraph, for each of the J\f chunks we construct a compressed kernel (table) 
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which approximates (j)e{iy, pa, Pa+i) as a sum of poles ^f{iy,pa,Pa+i)- The last step is to 
generate profiles (for the compression algorithm) on a large y-grid associated with a wide 
approximation window and sufficient resolution near the origin. However, these evaluations 
are now done via combination of all A/" tables. Therefore, they are drastically faster, since 
they are carried out through auxiliary evaluations made with the TV distinct pole sums 
(rather than ODE integration). Finally, we note that the physical teleportation kernel used 
in a numerical simulation is 

Sf (t,ri,r2) = (V(2M),ri/(2M),r2/(2M)), (57) 

where ^f{T,pi,p2) is the inverse Laplace transform (with respect to a) of ^(a, pi,p2)- 



4- Efficiency and storage 



As either £ or e~ becomes large the scaling observed in [41[ for compressed, blackhole, 
RBC kernels appears similar to the flatspace result (138|) described in Sec. IIII A4[ However, 
we stress that these are empirical observations, and there is no corresponding proof of 
( 138|) or a similar result in the blackhole case. Nevertheless, provided that the number of 
approximating exponentials in the blackhole case indeed grows sublinearly with both i and 
the inverse 1/e of the relative approximation error [cf. Eq. (149|) ]. our implementation of RBC 
satisfies the same efficient scalings (15^ estabhshed for flatspace compressed kernels. 

One might similarly ask whether or not our implementation of AWE for blackholes sat- 
isfles these scalings; we do not have an answer for this question, but here consider kernels 
from the i = 64 teleportation experiment considered later in Section IIV A[ Let us focus 
on the compressed (frequency domain) kernels ^eii^cr, 15) and ^64 (<7) 15, 240), respectively for 
RBC at pb = 15 and teleportation from pi = 15 to p2 = 240. Figure S] has already depicted 
the proflles from which the 25-pole approximation ^Q4{a, 15) is constructed. Notice that 
the approximation ^64 (c", 15, 240) has 32 poles, whereas the corresponding exact flatspace 
teleportation kernel would have 64 poles. Similar reduction for the £ = 64 occurred for 
compressed flatspace kernels considered earlier. As depicted in Fig. [8] and similar to the 
situation encountered in Fig. [3l the pole locations for ^Q^^a, 15) and ^64(c"; 15,240) are dif- 
ferent. Finally, we remark that we have tried to enforce the condition that the teleportation 
kernel has the same pole locations as the RBC kernel [cf. the discussion around Eq. Q]; 
however, we are then unable to achieve a compressed kernel with any accuracy whatsoever. 

The previous paragraph has considered the large-£, small-e limits. However, in this 
paper we mostly consider e fixed (typically machine precision) and small i, in which case, 
as we have seen, d > i. We remark that this situation is similar to the case of ordinary 
wave propagation on 2 + 1 flat spacetime. In that setting the \ow-n (Fourier index) "circle 
kernels" are also expensive to evaluate, with scalings similar to (138|1 only exhibited in the 



large n limit [49 



While the previous discussion has presented scalings for various limits, in practice our 
implementation of RBC/teleportation for £ = 2,3, and 64 amounts to adding on the order 
of 20 to 30 points to the spatial domain, and modest increase in work and storage costs for 
the numerical experiments considered here. 



23 



X 1^64(0-, 15)1 poles 



o|e|4(a, 15, 240)1 poles! 
6 o 



10 



10 



"""l^64(iy,15)| 
-1^^4(1;/, 15, 240)1 



i 



-10 



Q O o 



-5 

RecT 




FIG. 8. Compressed RBC and teleportation kernels for (. = 64. The right pane plots the 
modulus of the kernels for positive y. For the RBC kernel (dotted line) this corresponds to the 
modulus of the combined profiles shown in Fig. HI but here with only half of the domain for the 
independent variable y. 



IV. NUMERICAL EXPERIMENTS 

To carry out numerical simulations, we have used both the nodal Legendre discontinuous 
Galerkin method described in Ref. (further details of this method will not be given here) 
and a nodal Chebyshev method. Both methods feature multiple subdomains and upwinding. 



A. Pulse teleportation 

First consider the ^ = 2 Regge- Wheeler equation with the same initial data ( !T0|) given 
in Subsection III CI Using our multidomain nodal Chebyshev method, we perform five sep- 
arate evolutions on domains with outer boundaries taken as the h values corresponding to 
n = 30M, 60M, 120M, 240M, and 480M. We have respectively used 32, 37, 45, 62, and 
95 subintervals of uniform size, and in each case with 32 Chebyshev-Lobatto points per 
subinterval. Therefore, the spatial resolution for each evolution is comparable to the oth- 
ers. Evolutions are performed by the classical 4-stage explicit Runge Kutta method with 
timestep At ~ M(2.6794 x 10~^). For each evolution the inner boundary is a = — 200M, 
and therefore the Sommerfeld boundary condition — (H + $) = at the inner boundary 
is essentially exact. For all choices of outer boundary b we adopt the Laplace convolution 
RBC dH). Tables for = 30M, 60M, 120M, 240M, and 480M respectively have 19, 19, 19, 
18, and 17 poles, with each table computed in quadruple precision to satisfy the tolerance 
e = 10"^^. These tables are available at [3]. 

In all cases the field \E'(t,rb) is recorded as a time series at the boundary b, and in all 
cases but the last [b corresponding to r?, = 480M) we "teleport" the field from ri = 
to r2 = 480M. Each approximate teleportation kernel Sf ~ $2 features the same pole 



24 



10" 



10" 



10" 



10" 



1 1 1 1 1 


1 1 1 


read-off at 30M - 

— read-off at 60M , 

— read-off at 120M ■ 

— read-off at 240M 
1 1 



500 



520 



500 



520 



540 



560 




tiM 



580 



— teleported 30M - 


^ 480M 


— teleported 60M - 


>480M ' 


— teleported 120M - 


480M 


— teleported 240M - 


480M 


540 560 


580 



FIG. 9. Errors in read-off and teleported i = 2 waveforms relative to read-off 

WAVEFORM AT r2 = 480M. 



locations as the corresponding approximate RBC kernel S ^ f22. For the last = 480M 
simulation we simply record the field at the boundary, with this record then serving as 
a reference time series. We account for time delays by starting all recorded times series 
(whether read off or teleported) at time h — 6M. 

The top panel in Figure plots the errors in the waveforms recorded at the different h 
boundaries as compared to the reference = 480M waveform; as expected the systematic 
errors are large. The bottom panel plots the errors in the (ri = to r2 = 480M) teleported 
time series relative to the reference time series. With the reference time series viewed as 
the "asymptotic signal" , this "AWE" clearly yields 10 or more digits of accuracy relative to 
simple read-off. We have found similar results using other "pulses" based on polynomial, 
Lorentzian, and trigonometric profiles (in all case with the initial data initially supported 
away from the boundary, either exactly or to machine precision). 

We repeat the experiment for two different £ values. First, for £ = 3 we adopt the same 
initial data and experimental setup, except for the numerical tables which specify the RBC 
and teleportation kernels. For the £ = 3 experiment the number of poles for the RBC 
tables is either 15 or 16, and the number of poles for teleportation tables ranges from 15 
to 20. The results, shown in Fig. [TOl are comparable to those for (i = 2. Lastly, we repeat 
the experiment for £ = 64. For such a high ^ the evolutions are much more expensive, 
due to finer oscillations in both space and time. We now use 42 points per subdomain, 
with the number of subdomains typically increased by a factor of 2 or 3 relative to the 
numbers given above for £ = 2. Moreover, we adopt the timestep At = M(4.0461 x 10~^) 
and inner boundary a = — 80M which is casually disconnected from each waveform read- 
off/teleportation at the outer boundary. For i = 64 our RBC tables have between 23 and 25 
poles, and our teleportation tables between 30 and 32 poles. While these tables are large, 
note that even the corresponding exact flatspace RBC kernels would have 64 poles. Hence, 
in this experiment the savings afforded by kernel compression is already evident. Results 
are depicted in Fig. [TTJ 
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64 WAVEFORMS RELATIVE TO READ- OFF 



B. Luminosities from extreme- mass- ratio binaries 



An extreme mass ratio binary (EMRB) is a system comprised of a small mass-m^ com- 
pact object (the "particle") orbiting a much larger mass-M blackhole, where the mass ratio 
rrip/M <^ 1. EMRB systems are expected to emit gravitational radiation in a low frequency 
band (10~^ to 10~^ Hz), and therefore offer the promise of detection by a space-based grav- 
itational wave observatory like the earlier proposed LISA project j58|. Located within the 
solar system, such an observatory would be well approximated as positioned at J^"*" relative 
to expected sources. 
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A standard method for studying EMRBs uses the perturbation theory of Schwarzschild 
blackholes in an approximation which treats the particle as a point-hke Dirac delta function. 
The particle follows a timelike geodesic in the background Schwarzschild spacetime and is 
responsible for generating small metric perturbations which radiate away (see Refs. [l7- lol 



for modern accounts of the subject). Here we note that the axial metric perturbations for 
each (£, m) mode may be combined to form a gauge invariant scalar quantity "^f^^ which 
obeys Eq. O with the Regge- Wheeler potential \/^^(r) and a distributional forcing term 
Sf^^{t, r)u Likewise, the polar metric perturbations for each {£, m) mode may be combined 
to form a gauge invariant scalar quantity \l/f^ which obeys Eq. ([1]) with the Zerilli potential 
V'^{r) and a distributional forcing term Sf^{t, r). Both Sf^ and Sf^^ are built from linear 
combinations of an [i, m) mode decomposition of the stress-energy tensor and, as such, 
depend on the particle trajectory {rp{t),Tc/2, (f)p{t)) in the equatorial plane. Bounded and 
stable orbits are characterized by an eccentricity constant e and a semi-latus rectum constant 
p. Upon specification of (e,p), the resul ting trajectory is found by integrating the relevant 



system of ODEs given in Eq. (5) of Ref. [40|. Appendix C of [40] gives exact expressions for 
Sf™"'^ (see also [Tol. [H?!). 

With both 5'^'^'^ specified, we numerically solve for ^I/^^'^, starting with trivial initial 
data and smoothly turning on the source term Sern{t,T) over a timescale r ~ 150M to 
450M to prohibit static Jost junk solutions which may appear in some formulations when 
using inconsistent initial data (goI . Respectively, Sommerfeld and Laplace convolution 



RBC are enforced at the left and right physical boundaries (cf. Sec. Ill Aj) . The computational 
domain is given by the interval [— 400M, b], where the tortoise coordinate value b corresponds 
to rj, = 60M in Schwarzschild radius. Notice that as an approximation to the asymptotic 
signal at J^~^, the waveform read-off at will have an 0(r^^) systematic error, suggesting 
relative errors greater than one percent for = 60M. For our simulations, we have chosen 
16 and 3 subdomains to the left and right of the delta function respectively and represent 
the numerical solution by an order-40 or order-46 polynomial on each subdomain. The 
distributional source terms Sf^^''^ determine jump conditions in the fields at the 



particle location which we impose as junction conditions between subintervals |40[, with the 
motion of the particle incorporated through a time-dependent coordinate transformation. 
Our particular choice of Vf, = 60M ensures that the particle does not come too close to 
the outer computational boundary which might lead to over stretching of the coordinates. 
Temporal integration is carried out with a classical fourth-order Runge-Kutta method with 
timestep At = M(5 x 10"^). 

Computation of the luminosities for a particular orbit of the perturbing particle is a stan- 
dard benchmark test. For each (£, m) mode the energy and angular momentum luminosities 
at J^"*" (denoted by oo) and at the event horizon r = 2M (denoted by H) are given by 

+ (58a) 



oo,H 



L 



1 {i 


+ 2)! 


647r {i 


-2)! 


im {i 


+ 2)! 


647r {i 


-2)! 



^ We use the Cunningham-Price-Moncrief (CPM) master function [59.] which yields formulas for the axial 
sector which are on the same footing as those for the polar sector. 
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where the average {A{t)) of a time series A{t) is computed as 

1 



(A) 



A{t)dt, 



T2 - Ti = T,. 



(59) 



Here is the period of radial oscillation for the particle orbit. 

Before presenting our numerical results, we remark on the potential sources of error in an 
evaluated asymptotic waveform. At b we record both ^'^^'^'^(t, and their first temporal 
derivatives as a time-series. With the numerical setup described above, the relative pointwise 
error associated with these read-off waveforms is better than 10~^°. An additional source 
of (systematic) error is due to trivial initial data which, both incorrect and inconsistent, is 
known to generate spurious junk. At a finite and fixed radial location, spurious junk radia- 
tion propagates away (the potential for static Jost junk is discussed in Ref. [s^), although 
due to backscattering a "junk error tail" may develop which decays more slowly. Tail fields 
are expected to fall off like at J^"*", and at a fixed (much smaller) radial value. Evi- 
dently, the situation is worse at J^~^ where junk error tails decay more slowly. Additionally, 
we often need to average luminosity quantities over long periods of time. Taken together, 
these facts conspire to make the temporal average of a J^^ waveform especially prone to 
contamination by junk error tails, even at late times and especially when high accuracy 
is desired. Unfortunately, simply waiting for junk errors to die out may not be practical, 
because ODE and PDE numerical integrators typically introduce numerical errors which 
grow linearly with time. While convolution with an approximate teleportation kernel will 
introduce additional error, we believe that the dominant errors in our asymptotic waveforms 
stem from numerical method error and spurious junk. 



m Alg. 






f 00 




FR 
AWE 


1.27486196317 xlO"* 
1.27486196187 xlO~* 


1.66171571270 xlO"* 
1.66171571269 xlO"* 










1 FR 
AWE 


1.15338054092 xlO^*^ 
1.15338054091 xlO"*^ 


3.08063328605 xlO"^ 

3.08063328606 xlO"^ 


1.44066000650 xlQ-^ 
1.44066000619 xlO"^ 


2.77518962557 xlO''^ 

2.77518962558 xlO~s 


2 FR 
AWE 


1.55967717209 xlO"* 
1.55967717211 xlO"* 


1.84497995136 xlO"^ 
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TABLE I. Mode- by- mode 1 = 2 luminosities for the eccentric orbit described in the text. For a 
particle of mass nip these values should be scaled by m^. The table compares our asymptotic- 
waveform evaluation (AWE) method with the accurate frequency domain (FR) luminosities. FR 
results refer to Table III of Ref. 571] and are quoted to a relative error of 10~^^. For this experiment 
the outer boundary is = 60M. 



Through — )■ Too = 2M(1 x 10^^) teleportation, we approximately obtain the signals 
^fen'^'^ at and with these signals compute the energy and angular momentum lumi- 
nosities (l58k .b). The orbital parameters (e,p) = (0.764124,8.75455) and initial location 
(rp,0p) = ((pM)/(l -|- e cos(7r/2)), 0) specify the particle's path. As described above, we 
slowly turn-on the distributional source over a timescale of r ~ 150M to 450M. A phys- 
ically meaningful luminosity measurement will not depend on our choice of r, and from 
this consideration we find that by 3000M the spurious junk's effect is minimal. Table |T] 
compares the i = 2 luminosity measurements at J^~^ with the accurate frequency-domain 
results reported in Ref. j57|. We match their stated accuracy to better than 9 digits. 

As our final experiment we consider a circular orbit specified by the orbital parameters 
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{e,p) = (0,10). Circular-orbit luminosity measurements are time-independent, thereby al- 
lowing us to (i) better understand the influence of junk error tails on J^"*" waveforms and 
(ii) estimate errors due to our AWE procedure in a clean setting. With the same numerical 
parameters used for our eccentric orbit simulations, we compute i = 2 luminosities at 
and compare them to accurate frequency- domain results generated with the code described 



in Ref. [57|. By T = 6000M our results agree with the frequency- domain results to within 
a relative difference of less than 10~^^. Furthermore, we flnd the same level of agreement 
when the outer boundary is moved inward to = 30M, in which case we use the tele- 
portation/RBC kernel tables given in Appendix [D] (with the longer teleportation table for 
AWE). 

As the flnal measurement time is taken earlier, the agreement becomes progressively 
worse due to spurious junk radiation. Indeed, the solid black line in Fig. [12] (left) plots the 



relative error |£'22(^) ~-E^22,frI/ -E22,fr a time-series, where E^^pii frequency-domain 
value. For comparison we also compute a "luminosity" quantity0 £'22('^) from ^fgl^;^)- The 
sohd black line in Fig. [12] (left) shows the relative error |i;^2(^) -^22 (6500M) 1/^22 (6500M), 

where E22(6500M) is a late-time computation less contaminated by spurious junk. Compar- 
ing the black and red lines, we see that the junk error tails at J^^ persist longer than those at 
the outer boundary b. This observation suggests that spurious junk radiation is a stubborn 
problem for high accuracy studies. The right panel of Fig. [T^] indicates that the energy lumi- 
nosity errors (due to spurious junk) respectively decay as and for a flxed radial value 
and J^~^. If we view ^exact as either ^22(65 00, b) or then these rates are consistent with 
the expected decay rate for field error tails ^E'junk tail and the relationship for a numerically 

computed energy luminosity E OC Inexact + ^junktailP - l^exactP + 2Re(^exact^junk tail) • 
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TABLE II. Mode- by-mode i = 2 luminosities for the circular orbit described in the text. For a 
particle of mass nip these values should be scaled by nip. The table compares our asymptotic- 
waveform evaluation (AWE) results with frequency domain (FR) results computed by the code 
described in Ref. [571]. We are grateful to S. Hopper for generating these previously unpublished 



FR luminosity values. 



V. CONCLUSIONS 

In the context of the Regge- Wheeler /Zerilli (describing blackhole perturbations) and 
ordinary (describing acoustic phenomena) wave equations we have developed a procedure 
for obtaining the asymptotic far-fleld signal from a time-series recorded at a flnite radial 
value located beyond the spatial compact support of the initial data. Furthermore, we have 

^ At finite radial values, especially ones this small, £^22 is certainly not the energy radiated by the system. 
However, this value is computable and, furthermore, is theoretically (although perhaps not numerically) 
constant for circular orbits. Therefore, our intention here is to quantify the effect of junk error tails on 
its computation. Approximation of E'^, perhaps by extrapolation, might rely on such measurements. 
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FIG. 12. Energy luminosity errors due to spurious junk radiation. Both panels show 
-^^,fr|/^^,fr (red line, r^) and |i;^2(0 " ^22(65 00M)|/i;^2(6500M) (black line, r,,) as 
a time-series. Comparing the black and red lines, we see that the energy luminosity error (due to 
spurious junk) is particularly persistent at J^"*". See the text for further explanation. 

viewed asymptotic waveform evaluation as a limiting case of signal teleportation between 
two finite radial values. For each of these wave equations our steps are to (i) write down the 
exact relationship for teleportation in the Laplace frequency domain, (ii) approximate this 
relationship along the inversion contour (of the inverse Laplace transform) by a sum of simple 
poles, and (iii) then represent, through inversion, the asymptotic signal as a convolution 
[cf. Eq. (IH])] of the solution with time-domain kernel [cf. Eq. (JTj)] comprised of damped 
exponentials. A similar recipe might be used both to impose boundary conditions for and 
evaluate asymptotic waveforms from perturbations of a Kerr blackhole. The Teukolsky 
equation describing such perturbations is, like the cases treated in this paper, separable in 
the frequency domain. 

Through pre-computed numerical tables specifying each exponential's strength and damp- 
ing rate, we have demonstrated that accurate asymptotic waveform evaluation through tele- 
portation can be easily implemented. Our simulations based on these numerical tables 
correctly exhibit as the asymptotic decay rate for ^ = 2 tails. We have also performed 
generic-orbit, extreme- mass-ratio, binary simulations. From the solution recorded as close 
as 30M and 60M, we compute far-field luminosities which agree with accurate frequency 
domain results to a relative error of better than 10~^ (10~^^ for circular orbits). Our studies 
indicate that spurious junk radiation is particularly problematic for computations, be- 
cause far-field luminosity errors (due to spurious junk) decay at a slow rate. These results 
have been achieved without a compactification scheme to include J^"*" in the computational 
domain. Instead, they have relied on a Laplace convolution ([8]) which, decoupled from a 
numerical evolution, can also be carried out as a post-processing step on an existing time 
series. Finally, we have demonstrated effective signal teleportation between two finite radial 
locations, for £ = 2,3, 64 and with relative errors ^ 10^^°. These demonstrations are a pow- 
erful test of our AWE/teleportation method's accuracy as well as a practical sanity check of 
its implementation. 

As discussed in Sec. IIIIB41 our Laplace-convolution RBC and AWE methods would seem 
efficient from work and storage standpoints. Lower I kernels appear similar to the case of 
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ordinary wave propagation on 2 + 1 flat spacetime. In that setting the low-n (Fourier index) 
"circle kernels" are also expensive to evaluate to account for tail-like phenomena. In any 
case, for the i = 2,3, 64 cases considered in this paper a kernel's overall computational cost 
is roughly equivalent to to adding 20 to 30 points to the spatial domain. Moreover, it is 
non-intrusive, requiring no grid stretching or supplemental coordinate transformations, and 
may be carried out at any radial value beyond the spatial compact support of the initial 
data and sources. Finally, in their frequency-domain form, our kernels might be used to 
implement radiation boundary conditions and AWE in frequency- domain codes. Here we 
envision that the kernels would flrst undergo a "Wick rotation" prior to use. 

While the results of this paper are encouraging, we believe that further careful study is 
merited. First, the task of computing teleportation/AWE tables is daunting for a number of 
reasons. The main one is cost. Here we refer to the offline cost in generating a table, not the 
cost incurred by the user implementing AWE with such a table. As discussed in Sec. IIII B 31 
generation of AWE kernels costs upwards of 10^^ floating point operations. Moreover, since 
the cost is offline, with the resulting numerical table then "good for all time", we believe 
the process should be carried out in quadruple precision in order to achieve e = 10~^^ 
error tolerances Jl This further adds to the cost. Due to the difficulties associated with 
the computation of AWE tables, we have not adequately isolated all sources of error in 
their construction. A systematic and optimized procedure for computing kernels would 
greatly reduce the offiine costs. One possibility is an application specific quadrature rule. 
So far, we have employed the familiar Gauss-Kronrod rule, which is designed for high-order 
integration of polynomials. We might instead design a quadrature rule which is exact for 
the corresponding flatspace kernels 63 . 

We plan to construct a family of RBC and AWE tables for general use: Regge- Wheeler 
and Zerilli tables likely for 2 < i < 64, boundary radii pb = 15, 30, 60, 120, and an evaluation 
radius poo = 1 x 10^^ (or p^o = oo with a semi-inflnite quadrature rule, see Sec. IIII B 3\} . All 



kernels used in this paper, as well as these others, will be available at [44 
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Appendix A: Error estimates 

Suppose that we have an "exact" kernel B{t) and an associated Laplace convolution 

{B*m){t)= f B{t-t')m{t')dt'. (Al) 

We then have the following result for the relative convolution error associated with using an 
approximate kernel A{t) in place of B{t): 

\A{s) — B{s)\ 

\\A*m-B* ^||l,(o,oo) < sup ' ^ \\B * ^I|l2{o,oo), (A2) 

provided that B(s) 7^ holds for all s G iM. If this condition fails, we have instead 

\\A*^-B* ^|U,(o,oo) < (27r)-i/2 1^^^) _ ^^^^1 . (A3) 

seiR 

Before discussing their consequences, let us verify ( ]A2I) and ( 1A3I) . 

A Laplace convolution may be viewed as a Fourier convolution, that is 

rt 



I 

Jo 



/oo 
B{t - t')m{t')dt\ (A4) 
-00 



if we adopt that viewpoint that B{t) and '^{t) are causal functions, i.e. B{t) = B{t)9{t) and 
\E'(t) = '^{t)9{t), where 6{t) is the Heaviside step function. With this viewpoint, the Fourier 
transform of '^{t), for example, is 

= ^= / exp(-iwt)*(t)rft, (A5) 
v2tt Jo 

with the following formal relationship holding between the Fourier and Laplace transforms^ 

^(w) = ^=§(ia;). (A6) 
V 27r 

To establish flA2p . we view A * and i? * as Fourier convolutions of casual functions 
(each convolution is again a causal function) , and with the Parseval and Fourier convolution 



^ Although its neglection does not spoil the final estimate, the factor of 1/V27r was neglected on page 4156 
of Ref. 42] (and on pages 23 and 24 of |arXiv : gr-qc/0401001| 73). 
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theorems find that 



\A*-^ - B * '^\\l2{o, 



oo) = WA*"^ - B 

= \\M - B^\\L,iR) 

A{uj) - B{uj) 



< sup 



sup 



Biu) 
A{s) - B{s) 



15^1 



L2{ 



(AT) 



Bis) 



L2( 



Using the inverse Fourier transform on the final term to work backwards (again with the 
Parseval and Fourier convolution theorems), we obtain (\A2\i . The alternative estimate (lASp 
follows by nearly the same calculation. 

We view flA2p as an estimate for either of the error quantities 



£i = ||n^(-,rfe) * '^e{-,rb) - Ei{-,rb) * r;,) || 2.2(0,00) 

S2 = \\^e{-,rh,r^) * ^'^(■,rb) - Sf(-,rb, Too) * ^Ei-,rb)\\L2{o,co), 



and (lASp as an estimate for 



(A8a) 
(A8b) 



(A9) 



Quantities (lASb .b) measure the quality of our numerical approximations for RBC and tele- 
portation kernels, while (1A9P measures the quality of using the exact waveform at a large 
(but finite) radius Too 7^ 00 as an approximation for the asymptotic waveform at c^~^. We 
only present details for flA8b ) and flA9p . 

To use the estimate ([M]) for (lA8k). let B{t) = ^i{t,rb) and A{t) = Ee{t,rb). Using 
Algorithm m we approximate B{iy) = ^le{iy, r) along the axis of imaginary Laplace frequency, 
demanding that 

\hi{iy,r) -Ee{iy,r)\ ^ _ ^^^^^ 



sup 



|^£(i|/,r) 



< e. 



Here 2^(1?/, r) = Aiiy) is a sum of simple poles specified by one of our numerical tables, and 
e is the desired tolerance. This formula is essentially ( H9l) written earlier in dimensionless 
variables for the blackhole case. With the above identifications, Eqs. (lAlOp and flA2p yield 



Ex < e\\Qe{-,rb) * ^^(■, rfe)||L2(o,oo)- 



(All) 



Note that (lAlOp is an a posteriori bound; it is verified in step 5 of Algorithm |H 

To facilitate the analysis for (lAQp . we use ri and r2 in place of rb and r^o, with B(t) = 
$£(t,ri,oo) and A(t) = $£(t,ri,r2). Now referring to the absolute estimate (lASp . we must 
control the factor 



\Ais)-Bis) 



W{sr2,2Ms) 



W{sri,2Ms) W{sri,2Ms) 



(A12) 



here expressed for the blackhole case [cf. Eq. f l53|) ]. Further analysis of the blackhole 
case would presumably rely on the known asymptotic expansions (see footnote HJ for 
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W{sr,2Ms) = W{ap,a), but would seem difficult. Therefore, to proceed, we switch to the 
simpler flatspace case by setting M = and using the formal result W{sr, 0) = W{sr). The 
expression (fT5|) for IV^(sr) then gives 



\Ais) 

We now show that 



B(s) 



-k 



t-k 



ri 

sup I — 

sGiR V^2 



Y!k=i(^f^k{sr2 



\l-k 



0(l/r2). 



(A13) 



(A14) 



To establish this claim, note that the denominator of the expression inside the operation 
of complex modulus is the Bessel polynomial with zeros hi^k- Therefore, we expand the 
expression as a sum of simple poles, thereby finding 



sup 



Y!k=iCik{sr2 



< 



sup^ 



f^jiri,r2) 



sri - bej 



f^jiri,r2) 



Rebf 



(A15) 



The residue formula for a simple pole shows that each coefficient in the expansion obeys 
fij{ri,r2) = 0{l/r2), establishing the claim. These calculations show that (returning to Vb 
and Too for ri and 



^3 = 0(l/roo)||^£(-,r,) 



IL2(0,oo)" 



(A16) 



While we have not proved a similar formula for the blackhole case, this formula (an a priori 
estimate) has motivated our choice r^o = 2M(1 x 10^^) for double precision arithmetic. 



Appendix B: Derivation of NRBC without Laplace transform 



This appendix derives the nonlocal nonreflecting boundary condition (12^ without ap- 
pealing to the Laplace transform. In order to elucidate the main ideas, we choose to focus 
on the representative i = 2 case; the derivation for generic i is more cumbersome but similar. 
The Sommerfeld residual for the solution (IT^ is 



^^^2 + dr^2 = -^f^'\t - r) - ^f{t - r). 



(Bl) 



and we will show that this equation can be expressed as 

dt^2 + dr^2 = - / ^2it- t', r)^2(t', r)dt' (B2) 
r Jo 

+ -s(^)[4/«(-r) + 4/(-r)j}, 
\ Zr / ir r j J 



--(-|){-(f)[^/'"<-") 

where the i = 2 time-domain kernel is [cf. ([23 



\l2\t1 rj = — exp — t H exp — t 

r \ r / r \ r 



3 ^.\/3 



(B3) 
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Let us postpone establishing f lB2p . and first consider its consequences. If we assume ([H]) 
and evaluate (lB2p at r = r^, then the last two terms on the righthand side vanish and we 
have the desired result 

[dt^2it, r) + 9,vl>2(t, r)] l^^^ = 1 [n^i; n) * vl>2(., n)] (t). (B4) 
Indeed, from (i) the identity 

fit -r) = Ir^dt + dryidt + dr)^2{t, r) (B5) 



and the assumptions that (ii) \E'2(0,r) = = (9t\E'2)(0, r) in a neighborhood of and (iii) 
\1'2 obeys the radial wave equation, we conclude that f^^\—rb) = = f{—rb). Therefore, as 
claimed, the i = 2 case of holds exactly subject to our assumption f|T4l) on the initial 
data. Note that, even if f|T^ does not hold, the last term in flB2p decays exponentially. 



Now let us verify ( IB2p . assuming only the outgoing solution f[T^ without any restriction 

on the initial data. That is, we only assume that \l/2 = ^1^2^'*, with no contribution from \E'2~^'* 
[cf. ( ITT]) ]. First consider the quadratic polynomial 



z^ + 3z + 3, p{z±) = 0, 



(B6) 



where K^{z) is the MacDonald function (cf. footnote [2]). We appeal to the form (fT3l) of \l/2, 
and via integration by parts shift all time derivatives off of / and onto the exponentials. For 
generic z (either z^ or Z-) the calculation gives 



exp 



-it-t') 



f'\t'-r) + -f'\t'~r) + -J{t'-r) 



dt' 



— exp ( -t 



fW(t-r) + ^^{z' + 3z)f{t-r) 



^f^'\-r) + -,{z' + Sz)fi-r) 

■t 



+ (^2 + 3^ + 3) 



z 



exp 



z 



Lr 



{t-t') fit'-r)dt'. 



(B7) 



Since z is either z+ or Z-, the prefactor in the last term is p{z±) = 0. Therefore, 



exp 



-{t-t') ^2(t',r)dt' 



exp 



L r 



-{t-t') ^2{t',r)dt' 



-^{z^ + z_)f^^\t -r) + -(4 + 3z+ + zl + 3z_)f{t - r) 



r 

— exp ( —t 

\ r 

— exp ( — t 

\ r 



)fe/^^H-r) + ^(4 + 3z,)/(-r) 



(B8) 
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Combination of the last result with z± = — 1(3 =F i"\/3) then gives 
I f ^2{t - t', r)^2{t', r)dt' = -yj^'\t - r) - ^/(t - r) 

from which we immediately get (lB2p . 



(B9) 



( 



\ 2r 



Appendix C: RBC for other foliations of Schwarzschild 

In terms of the standard time slices and area radius the Schwarzschild line-element is 

ds^ = -fde + f-^dr^ + r'^de'^ + sin^ 0ci0^ f = 1- 2M/r. (CI) 
Define the outgoing (future and outward pointing) null vector 

^ ^ ^ + =^ f'^'^^ = 4 + (C2) 



/^^ dt dr dt dx ' 

again where x = is the tortoise coordinate. The exact RBC for these coordinates is 
essentially ( 145|) . with appropriate rescalings by 2M factors. In particular, with i7^(t, r) = 
(2M)-ia;^(t/(2M),r/(2M)), we write the RBC as 

/i/2^+[v^]=r-V(^^£*^), (C3) 

where = X from ([2]). To implement the boundary condition, we approximate it 

through the replacement — )■ S^. 

Following Zenginog lu's Eol analysis, we now consider a change of time slices defined by 
the new time variable 

A = t-/i(r), (C4) 
where hir) is the height junction. In terms of A the line-element becomes 

ds^ = -N^dX^ + Qrridr + VdX){dr + V'dX) + r^dO"^ + sin^ (C5) 

where the lapse, radial lapse, and radial component of the shift vector are respectively 



r2 f r- ^ 



^= l_lfHr ^=iV' ^' = -fHN\ (C6) 

Here H = dh/dr is the derivative of the height function. Define the outgoing (+) and 
incoming (— ) null vectors 
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Then z'^ = exp{-&)w^ and = exp{—i))z~^, where the boost angle is 



1 



log 



1 + ^N-'V 



log 



1 - y^A^-^l/'- 
Therefore, with respect to the new slices the exact RBC is 



1 - NfH^l^UW 
1 + NfH^l - UHf 



(C8) 



(C9) 



and it can similarly be approximated through the replacement VLi S^. As given by 
Zenginoglu [30|, the H functions for ingoing Eddington-Finkelstein and constant mean cur- 
vature foliations are respectively 



2M 



iEF 



r - 2M' 



J 



CMC 



(CIO) 



where J = ^Kr — Cr ^ in terms of the trace K of the extrinsic curvature tensor (based on 
Wald's definition 62| of the tensor) and constant C of integration. 



Appendix D: Numerical Tables 

This appendix collects the tables used for the numerical simulation documented in Sub- 
section [ircl Table UTTl determines the kernel S2(t, 30M) which approximates the exact ker- 
nel fi2(t,30M) = {2M)-^uj2it/{2M),15). The 19 locations /32,g and strengths 72,5 which 
make up this table have been computed in quadruple precision and satisfy the tolerance 
e = 10~^^. Entries of . OOe+00 correspond to outputs from the Alpert-Greengard-Hagstrom 
compression algorithm which are typically in the range 10~™ to 10"^*^°. 

We provide two different approximations for the time-domain teleportation kernel 
$2(t, 30M, 2M(1 X lO^^)) = (2M)-V2(t/(2M), 15, Ix 10^^), each denoted Sf (t, 30M, 2M(lx 
10^^)). Table [IV] determines the first Hf (t, 30M, 2M(1 x 10^^)). For this table notice that 
the 19 locations exactly match the /32,g listed in Table lllli Therefore, with this table 
the teleportation can be performed without evolving supplemental convolutions. However, 
we believe the tolerance for this table is only e = 5 x 10~^. Table IVl determines the second 
Sf (t,30M, 2M(1 X 10^^)) which now has 26 locations jS^g and strengths 7^^. Use of this 
table for teleportation with the RBC specified by Table HTll requires the evolution of 26 extra 
convolutions. However, we believe that this second approximate kernel satisfies a tolerance 
of e = 2 X 10-1^ 



[1] E. W. Leaver, Solutions to a generalized spheroidal wave equation: Teukolsky's equations 
in general relativity, and the two-center problem in molecular quantum mechanics, J. Math. 
Phys. 27 (1986) 1238-1265. 

[2] L. Barack, Late time dynamics of scalar perturbations outside black holes. I. A shell toy model, 
Phys. Rev. D 59 (1999) 044016 (20 pages). 
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Regge-Wheeler RBC table for ell = 2 and rho =15.0 

Gamma strengths Beta locations 



-2 


6076002831928367e- 


-08 


+0 


00e-^00 




-5 


41465293414875816- 


-01 


+0 


00e-^00 


-1 


7937477396220654e- 


-06 


+0 


OOe-i-00 




-4 


13109549893964766- 


-01 


+0 


00e-^00 


-3 


2816441859083765e- 


-05 


+0 


OOe-i-00 




-3 


19113384820765576- 


-01 


+0 


00e-^00 


-2 


8179763264971427e- 


-04 


+0 


OOe-i-00 




-2 


47112198718996596- 


-01 


+0 


00e-^00 


-1 


4509759948015657e- 


-03 


+0 


OOe-i-00 




-1 


91081637229234716- 


-01 


+0 


00e-^00 


-4 


4918693070976545e- 


-03 


+0 


OOe-i-00 




-1 


47496015587184506- 


-01 


+0 


00e-^00 


-5 


6790046261682662e- 


-03 


+0 


OOe-i-00 




-1 


13662999459085886- 


-01 


+0 


00e-^00 


-2 


0012016782502274e- 


-03 


+0 


OOe-i-00 




-8 


64769353811643416- 


-02 


+0 


00e-^00 


-2 


9649254206011509e- 


-04 


+0 


OOe-i-00 




-6 


45120651754510366- 


-02 


+0 


00e-^00 


-3 


2913867328382246e- 


-05 


+0 


OOe-i-00 




-4 


73323744420445576- 


-02 


+0 


00e-^00 


-3 


2675049152330702e- 


-06 


+0 


OOe-i-00 




-3 


41157754846636026- 


-02 


+0 


00e-^00 


-2 


8887585153331239e- 


-07 


+0 


00e^-00 




-2 


40489357047596546- 


-02 


+0 


00e-^00 


-2 


1640495893086479e- 


-08 


+0 


OOe-i-00 




-1 


64686329192834806- 


-02 


+0 


00e-^00 


-1 


2772861871474360e- 


-09 


+0 


OOe-i-00 




-1 


08456904230586966- 


-02 


+0 


00e-^00 


-5 


3164468909323526e- 


-11 


+0 


OOe-i-00 




-6 


75529185978649476- 


-03 


+0 


00e-^00 


-1 


2736896522814067e- 


-12 


+0 


00e-^00 




-3 


85256301968913256- 


-03 


+0 


00e-^00 


-1 


0598024220301938e- 


-14 


+0 


OOe-t-00 




-1 


84812150407888666- 


-03 


+0 


00e-^00 


-8 


9530431033189126e- 


-02 


+6 


2063746326002998e- 


-02 


-9 


47794908152390236- 


-02 


+5 


99279798774887206-02 


-8 


9530431033189126e- 


-02 


-6 


2063746326002998e- 


-02 


-9 


47794908152390236- 


-02 


-5 


99279798774887206-02 



TABLE III. Radiation boundary conditions. As indicated this table corresponds the ri, = 
30M and £ = 2. 



Regge-Wheeler extraction table for ell = 2 and rhol = 15.0 to rho2 = 1.0e-H5 
GammaE strengths BetaE locations 



-1 


75762630576795886- 


-08 


+0 


006-1-00 




-5 


41465293414875816- 


-01 


+0 


00e-^00 


-6 


41805142932012446- 


-08 


+0 


006-1-00 




-4 


13109549893964766- 


-01 


+0 


00e-^00 


-6 


27329710500936456- 


-06 


+0 


006-1-00 




-3 


19113384820765576- 


-01 


+0 


00e-^00 


-6 


93631179889879856- 


-05 


+0 


006-^00 




-2 


47112198718996596- 


-01 


+0 


00e-^00 


-5 


71806377507933456- 


-04 


+0 


006-^00 




-1 


91081637229234716- 


-01 


+0 


00e-^00 


-2 


78842475771758256- 


-03 


+0 


006-1-00 




-1 


47496015587184506- 


-01 


+0 


00e-^00 


-5 


88367920335704066- 


-03 


+0 


006-1-00 




-1 


13662999459085886- 


-01 


+0 


00e-^00 


-3 


65491361328921946- 


-03 


+0 


006-1-00 




-8 


64769353811643416- 


-02 


+0 


00e-^00 


-1 


04987467674996286- 


-03 


+0 


006-1-00 




-6 


45120651754510366- 


-02 


+0 


00e-^00 


-2 


42047818789951816- 


-04 


+0 


006-1-00 




-4 


73323744420445576- 


-02 


+0 


00e-^00 


-5 


57244641766299106- 


-05 


+0 


006-1-00 




-3 


41157754846636026- 


-02 


+0 


00e-^00 


-1 


21572967935489606- 


-05 


+0 


006-1-00 




-2 


40489357047596546- 


-02 


+0 


00e-^00 


-2 


66518132471934866- 


-06 


+0 


006-^00 




-1 


64686329192834806- 


-02 


+0 


00e-^00 


-4 


86617089811827696- 


-07 


+0 


006-1-00 




-1 


08456904230586966- 


-02 


+0 


00e-^00 


-8 


61836776120600446- 


-08 


+0 


006-^00 




-6 


75529185978649476- 


-03 


+0 


00e-^00 


-9 


37350711899108106- 


-09 


+0 


006-1-00 




-3 


85256301968913256- 


-03 


+0 


00e-^00 


-8 


78817870230940766- 


-10 


+0 


006-1-00 




-1 


84812150407888666- 


-03 


+0 


00e-^00 


-9 


11645360275914336- 


-02 


-5 


39537091551987806- 


-02 


-9 


47794908152390236- 


-02 


+5 


99279798774887206-02 


-9 


11645360275914336- 


-02 


^-5 


39537091551987806- 


-02 


-9 


47794908152390236- 


-02 


-5 


99279798774887206-02 



TABLE IV. Teleportation table for AWE. Note that the locations in this table match those 
in Table |TTI1 



[3] M. Piirrer, S. Husa, and P. C. Aichelburg, News from critical collapse: Bondi mass, tails, and 

quasinormal modes, Phys. Rev. D 71 (2005) 104005 (13 pages). 
[4] A. Zenginoglu, A hyperboloidal study of tail decay rates for scalar and Yang-Mills fields, 

Class. Quantum Grav. 25 (2008) 175013 (13 pages). 
[5] R. Sachs, Asymptotic Symmetries in Gravitational Theory, Phys. Rev. 128 (1962) 2851-2864. 
[6] R. Geroch in Asymptotic Structure of Space-Time, edited by F. P. Esposito and L. Witten 

(Plenum Press, New York, 1977). 
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Regge-Wheeler extraction table for ell = 2 and rhol = 15.0 to rho2 = l.Oe+15 
GajimiaE strengths BetaE locations 



-2 


7644898994070847e- 


-08 


-1-0 


00e-^00 




-4 


7566048766905883e- 


-01 


+0 


00e-^00 


-2 


0673535427889927e- 


-06 


+0 


OOe-i-00 




-3 


5108913199590891e- 


-01 


+0 


00e-^00 


-4 


2379338886728862e- 


-05 


+0 


OOe-i-00 




-2 


6424678470854152e- 


-01 


+0 


00e-^00 


-4 


3426207863682094e- 


-04 


+0 


OOe-i-00 




-2 


0004528999326085e- 


-01 


+0 


00e-^00 


-2 


5469795864949394e- 


-03 


+0 


OOe-i-00 




-1 


5187870767201261e- 


-01 


+0 


00e-^00 


-6 


0248402659888447e- 


-03 


+0 


OOe-i-00 




-1 


1566719641292256e- 


-01 


+0 


00e-^00 


-3 


8754717050848465e- 


-03 


+0 


OOe-i-00 




-8 


73646836524982210- 


-02 


+0 


00e-^00 


-1 


0839648026423402e- 


-03 


+0 


OOe-i-00 




-6 


49960457640380710- 


-02 


+0 


00e-^00 


-2 


5061924186508839e- 


-04 


+0 


OOe-i-00 




-4 


78863734859380640- 


-02 


+0 


00e-^00 


-5 


8503552767044972e- 


-05 


+0 


OOe-i-00 




-3 


4988356300604928O- 


-02 


+0 


00e-^00 


-1 


4014861186166239e- 


-05 


+0 


OOe-i-00 




-2 


5326730622776284O- 


-02 


+0 


00e-^00 


-3 


3862132702752687e- 


-06 


+0 


OOe-i-00 




-1 


8135330085424492O- 


-02 


+0 


00e-^00 


-8 


0926034459750484e- 


-07 


+0 


00e-^00 




-1 


2824248788413890O- 


-02 


+0 


00e-^00 


-1 


8800680938314195e- 


-07 


+0 


OOe-i-00 




-8 


9386729575681410O- 


-03 


+0 


00e-^00 


-4 


1796375696032426e- 


-08 


+0 


OOe-i-00 




-6 


1274292303408907O- 


-03 


+0 


00e-^00 


-8 


7558713932917318e- 


-09 


+0 


OOe-i-00 




-4 


11965477869428560- 


-03 


+0 


00e-^00 


-1 


7001775199292710e- 


-09 


+0 


OOe-i-00 




-2 


7071621335894381O- 


-03 


+0 


00e-^00 


-3 


0014940247349529e- 


-10 


+0 


OOe-i-00 




-1 


7308275437738592O- 


-03 


+0 


00e-^00 


-4 


7007801744202052e- 


-11 


+0 


OOe-i-00 




-1 


0699222812596474O- 


-03 


+0 


00e-^00 


-6 


3144518165051684e- 


-12 


+0 


00e-^00 




-6 


3367274917212608O- 


-04 


+0 


00e-^00 


-6 


9169662522379380e- 


-13 


+0 


OOe-i-00 




-3 


5455566302145149O- 


-04 


+0 


00e-^00 


-5 


6820400420323507e- 


-14 


+0 


00e-^00 




-1 


8296880288072223O- 


-04 


+0 


00e-^00 


-2 


9748099840520155e- 


-15 


+0 


OOe-i-00 




-8 


29998336321655450- 


-05 


+0 


00e-^00 


-6 


5083672189429277e- 


-17 


+0 


OOe-i-00 




-2 


9042514390090429O- 


-05 


+0 


00e-^00 


-9 


1164550073798437e- 


-02 


-5 


3953205902806563e- 


-02 


-9 


47794946592877550- 


-02 


+5 


9928005360963245© 


-9 


1164550073798437e- 


-02 


^-5 


3953205902806563e- 


-02 


-9 


47794946592877550- 


-02 


-5 


99280053609632450 



1-02 
1-02 



TABLE V. Teleportation table for AWE. This table is a more accurate approximation to 
the kernel also approximated by Table ITVl 
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